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ABSTRACT 

When verifying a sophisticated numerical code, it is a usual practice to compare the 
results with reliable solutions obtained by other means. This work provides such solu- 
tions for the wavelength dependent dust radiative transfer problem. We define a set of 
benchmark problems in spherical geometry and solve them by three radiative transfer 
codes which implement different numerical schemes. Results for the dust tempera- 
ture and emerging spectra agree to better than 0.1%, and can be used as benchmark 
solutions for the verification of the dust radiative transfer codes. 
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1 INTRODUCTION 

Dust is abundant in the universe and many astronomical 
objects are associated with it. Most notable examples are 
galaxies and stars, at almost all points in their evolution. 
Dust surrounding an embedded source scatters, absorbs, and 
reemits radiation originating from the source. This process- 
ing usually results in an overall shift of spectral energy dis- 
tribution to long, infrared (IR) wavelengths^. Embedded 
sources may be entirely obscured by the dust at optical 
wavelengths, and the only available information is obtained 
at IR wavelengths, invisible from earth until recently. 

Thanks to the recent developments in IR techniques 
and facilities, the spectral energy distribution for many ob- 
jects is now available, and the amount and quality of data 
are steadily increasing. Interpretation of these data can offer 
insight to the nature of optically obscured objects, making 
their IR signature a powerful tool of analysis. However, be- 
cause of the complexity of the radiative transfer problem 
such analysis must be aided by sophisticated computational 
tools. 

The wavelength dependent dust radiative transfer prob- 
lem can be solved only numerically and early attempts to 
obtain approximate analytical and semi-analytical solutions 
assumed gray opacity (e.g. Chandrasekhar 1934, Kozirev 
1934). These early works were based on the Eddington ap- 
proximation that the radiation field is isotropic. The first 
direct numerical solutions were obtained by Hummer & Ry- 
bicki (1971) who employed iterations over variable Edding- 



* This is strictly true only in spherical geometry. Sources 
with toroidal, or disk-like dust distribution can show significant 
amounts of short wavelength radiation if viewed face-on. 



ton factor (the ratio of the second to the zeroth moment of 
the radiation intensity). However, they also assumed a gray 
opacity, an approximation too crude for detailed analysis of 
observations. 

The first calculations of the wavelength dependent 
dust radiative transfer in spherical geometry were per- 
formed by Scoville & Kwan (1976), and Leung (1976). 
Both attempts included some unrealistic assumptions and 
the first full approximation-free solution was obtained by 
Rowan- Robinson (1980, hereafter RR). His method is based 
on the direct integration of the radiative transfer equa- 
tion, also known as ray tracing. Yorke (1980) developed a 
method based on iterations over variable wavelength de- 
pendent Eddington factors, a non-gray extension of the 
method originally used by Hummer and Rybicki. Both ap- 
proaches remained the most widely used methods in subse- 
quent developments of new codes (e.g., Groenewegen 1993) 
whose number has been steadily increasing during last two 
decades. Gradually, the codes' capabilities have been ex- 
tended to two-dimensional geometries (e.g., Efstathiou & 
Rowan-Robinson 1990, Collison & Fix 1991, Men'shchikov 
& Henning 1997), and even to the three-dimensional case 
(Steinacker & Henning 1996). Yet, an important shortcom- 
ing remains in this development. 

When verifying a sophisticated numerical code, it is a 
usual practice to compare the results with an analytical so- 
lution, or with a reliable solution obtained by other means. 
There is no such solution available for the wavelength de- 
pendent continuum radiative transfer problem. Definition 
of such a radiative transfer problem involves several compli- 
cated functions as its input (e.g., the spectrum of the embed- 
ded source, wavelength dependent grain optical properties), 
making analytical solutions impossible. The only practical 
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approach is to compare solutions obtained by several in- 
dependent codes for a set of well defined problems. If all 
such codes produce the same results, then it is likely that 
these results are reliable, and can be used as solutions of the 
benchmark problems for code verification. 

This work defines a set of such benchmark problems 
in spherical geometry and compares solutions obtained by 
three radiative transfer codes which implement different nu- 
merical schemes. The benchmark problems are defined in 
Section 2, where we also briefly describe the radiative trans- 
fer codes which we used. Solutions for the dust temperature 
and emerging spectra are presented in Section 3, as well as 
the discussion of our results. 



2 BENCHMARK PROBLEMS 

Most of works describing infrared emission from astronomi- 
cal sources present various computational results. Yet, it is 
difficult to use those results for the code verification. While 
many of the model parameters are usually explicitly listed, 
very often it is not easy to recognize all the underlying as- 
sumptions, nor to reproduce the grain optical properties 
used in the calculations. To avoid these problems, all in- 
put properties for the benchmark problems presented here 
are defined analytically. Such approach will ease comparison 
of the results obtained by other codes with those presented 
here. 

2.1 Scaling Properties of the Radiative Transfer 
Problem 

Traditionally, detailed modeling of IR radiation involved nu- 
merous input quantities, necessitating a large number of cal- 
culations to obtain successful fits, and diluting the value of 
the resulting success. However, as pointed out by RR, the 
number of relevant parameters can be drastically reduced 
by employing the scaling properties of the radiative transfer 
problem. For example, the luminosity of the central source 
is irrelevant, a fact by and large ignored in modeling astro- 
nomical sources. 

The importance of scaling was recently emphasized by 
Ivezic & Elitzur (1995) who demonstrated that for given 
dust type, the IR emission from late-type stars can be suc- 
cessfully described by a single parameter — the overall op- 
tical depth -rfj All other quantities (luminosity, mass, mass- 
loss rate, etc.) enter only indirectly through their effect in 
determining r, and thus are irrelevant in the modeling of IR 
emission. It was subsequently recognized that this powerful 
scaling is a general property of radiative transfer and that 
it can be extended to arbitrary geometries (Ivezic & Elitzur 
1997, hereafter IE97). IE97 point out that the spectral shape 
is the only relevant property of the heating radiation when 
the inner boundary of the dusty region is controlled by dust 
sublimation. Similarly, the absolute scales of densities and 
distances are irrelevant; the geometry enters only through 
angles, relative thicknesses and aspect ratios. The actual 
magnitudes of densities and distances enter only through 
one independent parameter, the overall optical depth. Dust 

t Assuming a steady-state radiatively driven wind. 



properties enter only through dimensionless, normalized dis- 
tributions that describe the spatial variation of density and 
the wavelength dependence of scattering and absorption effi- 
ciencies. We now proceed to define the benchmark problems 
in terms of these fully scaled quantities. 

2.2 Definition of the Benchmark Problems 

A central point source is embedded in a spherically sym- 
metric dusty envelope with an inner cavity free of dust. The 
source radiates as a black body at a given temperature T*. 
The dust is in radiative equilibrium with the local radiation 
field, and this uniquely determines the dust temperature, Td, 
at every point in the envelope. The scale of Td is determined 
by Ti, the dust temperature at the inner boundary, r\. All 
radial positions can be scaled by this radius defining a new, 
dimensionless variable y = r/n. The dimensionless outer 
radius of the envelope, Y = n/ri, is a free parameter. The 
dust density variation with y is assumed to be a power law 
oc y~ p . The actual dust density, envelope size and opacity 
scales are all combined into a value for overall optical depth 
specified at some fiducial wavelength. We choose 1 fj.m for 
this wavelength and denote the corresponding optical depth 
by n. 

Dust optical properties are specified as scale-free wave- 
length dependent absorption and scattering opacities nor- 
malized to unity at 1 ^m, g abs = «A,abs/«i,ab8 and g sca = 
^A,sca/«:i,sca, respectively. For spherical amorphous grains 
with radius a, g a bs and g sca are roughly constant for A < 2na, 
and fall off as A -1 and A -4 for A > 2?Ya, respectively. To 
mimic this behavior^] with analytical functions, we choose 

gabs = <7sca = 1 (1) 

for A < l^rn, and 

gabs = J (2) 

gsca = (3) 

for A > lu.m. The arbitrary choice of 1 \im as transitional 
wavelength is motivated both by typical astronomical grain 
sizes, and by desire to have such a transition at short wave- 
lengths where its effects are most easily discernible. The cho- 
sen forms correspond to the maximal theoretically possible 
efficiencies for spherical amorphous grains with radius 0.16 
u.m (Greenberg 1971). 

The above listed quantities fully specify the benchmark 
problems. We perform calculations for T*= 2500 K, Ti= 
800 K and Y — 1000, and produce eight different models: 
for two different density distributions, p = 2 and p — 0, and 
four values of optical depth, n = 1, 10, 100, 1000. Selec- 
tion of different density distributions and a range in optical 
depth minimizes the chance that eventual inconsistencies in 
different numerical schemes would not be noticed. 

2.3 Code Description 

The three codes used in this work have similar approach to 
the solution of the radiative transfer problem, and none of 

T We assume isotropic scattering. 
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them introduces any approximations (e.g., closure relation in 
the moment methods, Auer 1984). Since the problem cannot 
be solved analytically, the solutions are obtained on discrete 
spatial and wavelength grids. Within arbitrary numerical 
accuracy which determines the grid sizes, the solutions can 
be considered exact. 

The spatial grids include the radial position and impact 
parameter grids, which also directly determine the angular 
grid needed for the integrations over solid angle. Sizes of 
the spatial grids range from ~10 to several hundred points, 
depending on the method and overall optical depth for a 
given model. The wavelength grid has typically ~100 points. 
To solve the radiative transfer problem means to determine 
the radiation intensity at each of the grid points, that is, to 
determine the wavelength dependent intensity as a function 
of angle for every radial position. 

The radiative transfer equation is an integro-differential 
equation; in other words, the intensity at any grid point 
depends on the intensities at all other grid points. Thus 
the problem cannot be solved by straightforward techniques, 
and various numerical methods rely on iterations of differ- 
ent types. These iterations can be performed either for the 
intensity, or for its moments (energy density, flux, pressure, 
etc.). The codes used in this work implement different nu- 
merical schemes to perform the iterations, and we proceed 
with their brief descriptions. 



Table 1. Dimcnsionlcss parameter ^ for eight benchmark models 



pO) n = 1 ri = 10 n = 100 n = 1000 

2 3.48 5.42 13.1 84.1 

2.99 3.00 3.10 3.75 



(a) Power for the power-law describing the dust density distribu- 
tion 



Code 1 This code was originally developed by Yorke 
(1980) and generalized by Men'shchikov & Henning (1997) 
and Szczerba et al. (1996, 1997). The original version solves 
differential moment equations obtained by analytically in- 
tegrating the radiative transfer equation multiplied by the 
powers of the direction cosine, over the solid angle (e.g. Auer 
1984). Such an approach always results in a number of dif- 
ferential equations of one less than the number of unknown 
quantites (the moments of the radiation intensity) . The clo- 
sure relation can be expressed in terms of the variable Ed- 
dington factor, and this code calculates it directly from its 
defining equation after every iteration. Later versions of the 
code were extended to improve explicit ray tracing, with it- 
erations repeated until convergence in dust temperature and 
mean intensity is achieved. 

Code 2 Groenewegen (1993) has developed a code 
which solves the radiative transfer equation in spherical ge- 
ometry from first principles, assuming isotropic scattering. 
The intensity is evaluated explicitly on sufficiently fine grids 
to obtain desired accuracy, and iterations are repeated un- 
til the convergence is achieved. This code was developed to 
allow for an explicit mass-loss rate dependence on time and 
velocity law, rather than a power-law density distribution. 
The model was tested against the results obtained with the 
model of Rogers & Martin (1984, 1986), for silicate grains 
up to an optical depth at 9.5 /im of 50; differences were 1% 
at most. For the purpose of the present paper, models with 
p — 2 (constant expansion velocity and mass loss rate) and 
ti = 1, 10, 100 were calculated (ri = 1000 was not included 
for this code because of unsatisfactory convergence). 

Code 3 The third code, DUSTY, was developed by 
Ivezic, Nenkova & Elitzur (1997) and is publicly available, ft 
solves the integral equation for the energy density obtained 



Figure 1. Dust temperature distribution through the envelope 
for two density distributions and optical depths at 1 um as 
marked. 

by analytically integrating the radiative transfer equation. 
The subsequent numerical integration is transformed into 
multiplication with a matrix of weight factors determined 
purely by the geometry. The energy density at every point 
is then determined by matrix inversion, obviating the need 
to iterate over the energy density itself. That is, unlike other 
codes, DUSTY can directly solve the pure scattering prob- 
lem. The intensity is not explicitly evaluated, and can be 
easily recovered from the source function by using the weight 
factor matrix. 



3 RESULTS 

The full solution of the radiative transfer problem is con- 
tained in the radiation intensity. However, even in the spher- 
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tained by different codes are smaiier than the fine thickness 
(< 0.1%). 

The detailed numerical values for emerging spectra pre- 
sented in Figures 1 and 2, and for all other relevant quan- 
tities can be obtained in computer readable form from Z. 
Ivezicpj. These results can be used for the verification of the 
wavelength dependent radiative transfer codes, especially 
for the new multi-dimensional ones. Also, as the numeri- 
cal radiative-hydrodynamical codes for modeling star and 
planet formation begin to include the wavelength depen- 
dent radiation transfer, the benchmark problems presented 
here might prove valuable for establishing confidence in the 
accuracy of their results (Boss 1996). 
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Figure 2. Spectral energy distribution of the emerging radia- 
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marked. 
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4<t7T 
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